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Abstract 

We consider a harmonically bound Brownian particle coupled to two distinct heat reservoirs at 
different temperatures. We show that the presence of a harmonic trap does not change the large 
deviation function from the case of a free Brownian particle discussed by Derrida and Brunet and 
Visco. Likewise, the Gallavotti-Cohen fluctuation theorem related to the entropy production at 
the heat sources remains in force. We support the analytical results with numerical simulations. 
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I. INTRODUCTION 



There is a strong current interest in the thermodynamics and statistical mechanics of small 
fluctuating non-equilibrium systems. The current focus stems from the recent possibility of 
direct manipulation of nano-systems and bio-molecules. These techniques permit direct 
experimental access to the probability distribution for the work and indirectly the heat dis- 
tribution [l-9]. These methods have also opened the way to the experimental verification of 
the recent fluctuation theorems, which relate the probability of observing entropy-generated 
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429] - 



191 ] which establishes a 



trajectories, with that of observing entropy-consuming trajectories 
We shall here focus on the Gallavotti- Cohen fluctuation theorem 
simple symmetry for the large deviation function \x for systems arbitrarily far from thermal 
equilibrium. Close to equilibrium linear response theory applies and the fluctuation theo- 
rem become s equiva lent to the usual fluctuation-dissipation theorem relating response and 
fluctuations 
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30]. 



A simple example of non-equilibrium system has been introduced recently by Derrida 
and Brunet [3JJ. In this model a particle or rod is coupled to two heat reservoirs at different 
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331 ] have shown that 



temperatures. We also note that Van den Broeck and co-workers 
an asymmetric object coupled to two heat reservoirs is able to rectify the random thermal 
fluctuations and thus exhibits a net motion along a preferred direction. It is therefore of 
interest to know whether the global behavior of these fluctuations, e.g., their fundamental 
symmetries, are left unaltered in the case one includes a potential or a particular interaction 
in such simple models. Furthermore, one is interested in knowing what type of interaction 
or lattice potential may increase, for example, the efficiency of a Brownian motor. When 
dealing with systems coupled to different heat baths, e.g., a chain of coupled oscillators, one 
of the main trends is to understand which essential properties of the microscopic dynamics 
lead to a diffusive limit for the energy 34|. Finally, it is also of importance to understand 
how heat conduction is affected when one deals with very small systems. 

More precisely, for a system driven into a steady non-equilibrium state by the coupling 
to for example two distinct heat reservoirs or thermostats at temperatures Ti and T 2 , a heat 
flux dQ/dt is generated in order to balance the energy. The heat flux is fluctuating and 
typically its mean value d(Q)/dt is proportional to the temperature difference. Focusing 
on the integrated heat flux, i.e., the heat Q(t) = J" * dr(dQ(T) jdr) over a time span t, 
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this quantity also fluctuates and typically grows linearly in time at large times. For the 
probability distribution we obtain the asymptotic long time behavior 



P{Q,t) oc e 



tF(Q/t) 



defining the large deviation function F(q). The Gallavotti-Cohen fluctuation theorem then 
establishes the symmetry 



F(q)-F(-q)=q[l/T 1 -l/T 2 



[1.2) 



Likewise, for the characteristic function 



(e AQ W) oc e^ x \ 



the fluctuation theorem states the symmetry relation 



;i.3) 



/x(A)=M-A + l/T 1 -l/T 2 



;i.4) 



The fluctuation theorem has been demonstrated under quite general and somewhat ab- 
stract conditions [l^ . It is therefore of importance to discuss the theorem in the context of 
specific models where the large deviation function //(A) can be derived explicitly. 

The large deviation function /z(A) can be determined explicitly for the simple non- 



equilibrium model introduced by Derrida and Brunet [31]; this model has also been discussed 
by Visco 35) and Farago 36[. The model consists of a single Brownian particle or rod cou- 
pled to two heat reservoirs at temperatures T\ and T 2 with associated damping constant T\ 
and T 2 . Here the heat Q is transported from one reservoir to the other via a single particle. 
These authors find that the large deviation function has the explicit form 



MA) = - 



r x + r 2 - jy\ + v 2 2 + 2r 1 r 2 (i - 2\t 1 + 2AT 2 - 2A 2 t\t 2 ) 



;i.5) 



This expression for /x(A) is consistent with the boundary condition /i(0) = following from 
(jl.3p and in accordance with the fluctuation theorem (ll.4p . i.e., //(A) = A+l/Ti — 1/T 2 ). 
For T\ = T 2 the large deviation function /x(A) is symmetric, i.e., //(A) = //(—A). In this case 
the heat fluctuates between the two reservoirs and there is no net mean current. If we 
decouple one of the reservoirs by setting T 2 = (or r x = 0) the system is in equilibrium 
with a single reservoir and we have /x(A) = for all A. Finally, from (I1.3P we infer the mean 



value (the first cumulant) and the second cumulant 



t 1 1 "hi 2 

(Q 2 ) - (Q} 2 _ 2T 1 T 2 T 1 T 2 2T|T|(7W^ 

t r! + r 2 (rx + r 2 ) 3 




(1.6) 



(1.7) 



Here we extend the Derrida-Brunet model to a Brownian particle moving in a harmonic 
trap and analyze the large deviation function. The paper is organized in the following 
manner. In Sec. [Ill we set up the model with focus on the heat transfer Q(t) and the large 
deviation function /i(A). In Sec. II III we evaluate the first and second cumulants within a 
Langevin approach, comment of the Fokker-Planck approach but focus in particular on the 
Derrida-Brunet method. We derive the differential equation for the characteristic function 
and determine the large deviation function. In Sec. IIVI we support the analytical findings 
by a numerical simulation. Sec. |V] is devoted to a summary and a discussion. 

II. MODEL 

We consider a ID Brownian particle harmonically coupled to a substrate by a force con- 
stant n. This configuration also corresponds to a Brownian particle in a harmonic trap. The 
particle is, moreover, in thermal contact with two distinct heat reservoirs at temperatures 
Ti and T 2 . The heat transferred in time t from the two heat reservoirs is denoted Q\ and 
Q 2 , respectively. Finally, the corresponding damping constants are denoted Ti and T 2 , re- 
spectively. The configuration is depicted in Fig. [IJ Denoting the position of the particle 
by u and the momentum by p and assuming m = 1, a conventional stochastic Langevin 
description yields the equation of motion 



du 



(2.1) 




(T 1 + T 2 )p -KU + & + 6, 



(2.2) 



dt 



where the Gaussian white noises £1 and £ 2 are correlated according to 



(mm) 



2TiTi6(t), 



(2.3) 



<6(*)6(o)> 



2T 2 T 2 6(t), 



(2.4) 



<6(*)6(o)> 



0. 



(2.5) 
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FIG. 1. We depict a harmonically bound particle interacting with heat reservoirs at temperatures 
Ti and T 2 . The heat transferred to the particle is denoted Q\ and Q 2 , respectively. The particle 
is attached to a substrate with a harmonic spring with force constant n. 

The heat flux from the reservoir at temperature Ti, i.e., the rate of work done by the 
stochastic force — F±p + £1 on the particle, is given by 



correspondingly, the heat flux from the reservoir at temperature T 2 has the form 

dQ 2 



(2.6) 



dt 



-r 2 p 



(2.7) 



The equations fl2.lti2.7p define the problem and the issue is to determine the asymptotic long 
time distribution for the transferred heats Qi and Q 2 , 



Q n (t) = / dr(-T n p 2 (r) +p(r)f„(r)), n = 1,2. 



(2.8) 

At long times the heat distribution in terms of its characteristic functions is given by ( II. 3p . 
i.e., 



(e AQ " w ) oc e' M " (A) , n = l,2, 



(2.9) 



where the large deviation function /i n (X) is associated with Q n (t). 



Noting that since the total noise £ = £i+£ 2 is correlated according to (£(£)£(0)) = (2FiTi+ 
2T2T2)5(t) and invoking the fluctuation-dissipation theorem {30J we readily infer that the 
system is in fact in equilibrium with the effective temperature T = (TiTi + r 2 T 2 )/(r 1 + T 2 ). 
This argument also implies that the stationary distributions for u and p are given by the 
Boltzmann-Gibbs expressions Po(p) oc exp(— p 2 /2T) and Po(u) oc exp(— nu 2 /2T). The non- 
equilibrium features are obtained by splitting the effective heat reservoir at temperature T 
in two distinct heat reservoirs at temperatures T\ and T 2 and monitoring the heat transfer. 
From the equations of motion (12. ip and (12.21) we infer two characteristic inverse lifetime in 
the system given by Yi + T 2 and n l l 2 . In the following we assume that the system is in a 
stationary non-equilibrium state at times much larger than (Ti + T 2 ) _1 and k^ 1 ^ 2 and thus 
ignore initial conditions, i.e., the preparation of the system. The role of the initial condition 
on the distribution P(Q,t) is a more technical issue, see Visco (35]. 



III. ANALYSIS 

We wish to address the issue to what extent the presence of the spring represented by 
the term ku in the equation of motion (12.21) changes the large deviation function (ll.5p in the 
free case. In the case of an extended system coupled to heat reservoirs at the edges, e.g., an 
harmonic chain, the heat is transported deterministically across the system and the large 
deviation function will depend on the internal structure of the system, e.g., in the harmonic 
chain the spring constant. For vanishing coupling the edges in contact with the reservoirs 
are disconnected and the large deviation function must vanish. However, for a single particle 
in a harmonic well there is no internal structure or internal degrees of freedom and the case 
is special. 

In addition to numerical simulations three analytical approaches are available in investi- 
gating this issue: i) a Langevin equation method taking its starting point in the equations 
of motion (I2.m2.2p and determining the distribution of the composite quantity Q n on the 
basis of a Greens function solution and Wick's theorem, ii) an analysis based on the Fokker- 
Planck equation for the joint distribution P(u,p,Q,t), and iii) a direct approach suggested 
by Derrida and Brunet which directly aims at determining the long time behavior of the 
characteristic function (exp(AQ(t)), yielding the large deviation function. 
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A. Langevin approach 



Here we delve into the Langevin approach and discuss the evaluation of the first two 
cumulants of the distribution P(Q,t). 



1. The first cumulant - The mean value 



The linear equations of motion (I2.lti2.2p readily yield to analysis. In Laplace space, 
ining u(s) = J °° dtu(t) exp(-st), etc., we obtain the solution 

p (s) = G(s)(Us)+Us)), (3-1) 



where the Greens function G(s), broken up in normal mode contributions, has the form 



G(.) = J2!_ + (3.2) 
s — Si s — s 2 



Here the resonances are given by 



Sl = -i[r + f], (3.3) 
s 2 = -\[T-n (3-4) 

r = r 1 + r 2 , (3.5) 



f = VT 2 - 4k; (3.6) 



we note the relations S\ + s 2 = — T, SiS 2 = K, and Si — s 2 = —y/T 2 + 4k. 
The the amplitudes mi and m 2 have the form 

mi = — — — , (3.7) 
si - s 2 

m 2 = (3.8) 

S2 - Si 

note the sum rule mi + m 2 . For T 2 > 4k the system is overdamped; for T 2 < 4k the 
system exhibits a damped oscillatory behavior with frequency \/2k — T 2 . In time we infer 
the solution 

p(*)= / dT( mi e Sl{t -^ +m 2 e s ^^)(Ur) + Ur))- (3.9) 
Jo 

We note that in the limit k — > 0, S\ — > —T, s 2 — > 0, mi — >■ 1, and m 2 — >■ 0, the position 
u is decoupled from the momentum p and we recover the model proposed by Derrida and 



Brunet 



3l|. 



Expressing time integration as a matrix multiplication and introducing the short hand 
notation p = (G\ + G2X6 + £2), where G n (t,t') = m n exp(s n (t — t'))r](t — t'), n = 1,2, we 
obtain from (I2.6H2.7I) 

^ = -r B ((G! + G a )(fc + 6)) 2 + ^n(Gx + G 2 )(6 + 6). (3.10) 
at 

For the mean flux d(Q n )/dt we then have averaging over the noises £1 and £2 according to 
(12.3112.51) 

= -r n (2r 1 T 1 + 2r 2 T 2 )(G 1 + g 2 ) 2 + 2iyr n (G 1 (o) + g 2 (o)). (3.11) 

Inserting / G n (t - t') 2 dt' = -m 2 n /2s n , J G x {t - t')G 2 (t - t')dt' = -mira 2 /(si + s 2 ), and 
G n (0) = m„r](0), 77(0) = 1/2, and reducing the expression we obtain 

- 2r n (r 1 r 1 + r 2 r 2 ) -J- + + — — + r n T n ( mi + m 2 ). (3.12) 



dt ' 2s 2 Si + s 2 j 

By insertion of m 1; m 2 ,Si, and s 2 the dependence on the spring constant k cancels out and 
we obtain 

(Qi> r x r 2 



t r\ + r 2 
{Q2) r x r 2 



(T1-T2), (3.13) 
■(T2-TO, (3.14) 



t ri + r 2 

independent of k and in agreement with the free particle case (11.61) . The independence of 
the mean value shows that the heat transport is unaffected by the presence of the spring. 
This feature is a result of the absence of internal structure in the single particle case. 



2. The second cumulant 



The evaluation of the second cumulant is more lengthy, involving Wick's theorem 37] 
applied to four noise variables. Focussing on Q = Q\ we have in matrix form 

(Q 2 ) = [ dr [ dT'a-r^GGZ + ^GOi-r^G'G'Z + ^G'O), (3.15) 
Jo Jo 

where G = Gi + G 2 , £ = 6+6, £GG£ = / dt'dt"£{t')G(T, t')G{r, t")f(t"), and £G'G'£ = 
/ dt'dt"i{t')G{r' ',t')G(T' ',*")£(*")• Applying Wick's theorem to the product (££££) entering 
in (I3.15P we note that only the pairwise contractions between the r and r' factors in (I3.15P 
contribute to the cumulant (Q 2 ) — (Q) 2 ', the contractions within the r and r' terms factorize 



in (I3.15P and yield (Q) 2 - Inserting £ = £1 + £2, applying Wick's theorem in pairing the noise 
variables, and using (I2.3ti2.5p . we obtain 

(Q 2 ) - (Q) 2 = fdt' f dt"[L{t',t") + M{t',t") + N{t',t")], (3.16) 
Jo Jo 

L{t't") = STl^Tt + T 2 T 2 ) 2 J drdr'Git', r)G{t\ r')G{t\ r)G{t", r'), (3.17) 

M(t't") = A(Tpf + T 1 T 2 T 1 T 2 )S(t' - t") J drG{t', r) 2 , (3.18) 

N{t't") = ST^TlT 2 + Y 1 Y 2 T 1 T 2 )G{t\t") J drG{t' \r)G{t" ,r). (3.19) 

Finally, inserting G = G\ + G 2 , using G n (t, t') = m n exp(s n (t — t'))r)(t — £'), and performing 
the integrations over t', t", r, and r', the dependence on the spring constant again cancels 
out and we obtain the free particle result 

(Q 2 )-(Q) 2 _ 2r 1 r 2 T 1 r 2 2T 2 t 2 (t 1 -t 2 ) 2 

t " r 1 + r 2 + (r 1 + r 2 )3 ■ (3J0) 

The Langevin approach turns out to be too cumbersome for the present purposes and we 
shall not pursue it further but note that the results for the two lowest cumulants corroborate 
the suggestion that the large deviation function is independent of the spring. 



B. Fokker-Planck approach 



Although we shall eventually complete the analysis using the Derrida-Brunet method, we 
include for the benefit of the reader and for completion the Fokker-Planck approach and the 
issues arising in this context. It is here convenient to consider the Fokker-Planck equation 
for the joint distribution P(u,p,Q,t), Q = Q\. It has the form 



dt 



dp 2 



dp 



_d_ 

dQ 



dP 



(p 2 + T 1 )P + T 1 p 2 — + 2T 1 p 



dQ 



dP 
dp 



where {P,H} denotes the Poisson bracket 

dP dH dP dH 



{P,H} 



KU- 



dP dP 
^ du 



(3.21) 



(3.22) 



dp du du dp dp 
The heat distribution after having analyzed the Fokker-Planck equation is then given by 



P(Q,t)= dudpP(u,p,Q,t). 



(3.23) 



Defining the characteristic function with respect to the heat by 

C(A) = J dQP(u,p,Q,t)e XQ , (3.24) 

and noting that d/dQ — > —A and d 2 /dQ 2 — > X 2 we obtain for C 

dC{\) 



dt 

where the Liouville operator L has the form 



L(A)C(A), (3.25) 



L{\)C{\) ={C(X),H} + (T x T a + T^)^! + ( ri + T 2 ) d{pC{X)) 



dp 2 dp 
(p 2 + T 1 )C(A) - XT lP 2 C(X) + 2T lP dC(X) 



(3.26) 



tip 

The case of an unbound particle Brownian particle for k = has been discussed in detail 



by Visco |35|, see also Farago [36|. Here {C(A), H} = —pdC(X)/du and integrating over the 
position u which is decoupled from the momentum p we obtain a second order differential 
equation for C of the Hermite type. By means of the transformation 

C (A) = e-W (5( A), ^(A) ^^; 2 ^' . (3.27) 

C(A) satisfies the Schrodinger equation for a harmonic oscillator and we infer the spectral 
representation 



C(A) = e -^)&> 2 -^) J2e En(X)t *n(p)*n(po), (3.28) 



n=0 



where —E n (X) is the discrete harmonic oscillator spectrum and ^ n {p) the associated normal- 
ized eigenfunctions. We have, moreover, imposed the initial condition C(t = 0) = S(p — po), 
where po is the initial momentum. The large deviation function is thus given by the ground 



state energy — E (X) yielding (jl.5P ; for further discussion see Visco [35 1. 

In the case of a bound Brownian particle for k ^ the Poisson bracket enters and the 
position of the particle comes into play. The Liouville operator becomes second order in u 
and p and is more difficult to analyze. We shall not pursue a further analysis of the Fokker- 
Planck equation here but anticipate, in view of the properties of the cumulants discussed 
above, that the maximal eigenvalue yielding \x remains independent of k. 
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C. Derrida-Brunet approach 



It is common to both the Langevin approach and the Fokker-Planck approach that they 
carry a large overhead in the sense that one addresses either the complete noise averaged solu- 
tion of the coupled equations of motion for u and p or the complete distribution P(u, p, Q, t). 



On the other hand, the method proposed by Derrida and Brunet [3l| circumvent these issues 
and directly addresses the large deviation function fi. 

Focussing again on Q = Q\ the long time structure of the heat characteristic function 

C(t) = (e AQW ) oce^ (A) , (3.29) 

immediately implies that C(t) satisfies the first order differential equation 

^ = MA)C(t). (3-30) 

The task is thus reduced to constructing this equation and in the process determine the 
large deviation function /i(A). 

In order to deal with the singular structure of the noise correlations as expressed in (12.31 



12. 5p and avoid issues related to stochastic differential equation 38j , it is convenient to coarse 
grain time on a scale given by the interval At and introduce coarse grained noise variables 

1 rt+At 

Fi = -^J Ur)dr, (3.31) 

I rt+At 

F2= AtJ t Ur)dT - (3 ' 32) 

Since £i and £2 are stationary random processes F\ and F 2 are time independent. Moreover, 
we have (Fx) = (F 2 ) = (FiF 2 ) = 0, and the correlations 

(F?) = (3.33) 
(Fi) = ™. (3.34) 

The coarse graining in time allows us to construct a difference equation for C(t) for then at 
the end letting At — > 0. Using the notation p(t + At) = p', etc., we thus obtain in coarse 
grained time from the equations of motion f l2.lti2.2p to O(At) 

u' = u + pAt, (3.35) 
p' =p+ (— (Tx + Y 2 )p -ku + F^ F 2 )At. (3.36) 

11 



For the heat increment we have from (12. 8p 



/t+At 
dr(-T 1 p(r) 2 +p(T)F 1 ), (3.37) 

Since from (I3.33P Fi is of order (At) -1 / 2 we must carry the expansion to 0((At 2 )) and we 
obtain 

Q' = Q + (F lP - T lP 2 )At + ^{F l F 2 + F 2 )(At) 2 . (3.38) 

We next proceed to derive a difference equation for C. This procedure will in general produce 
correlations of the type (e x ®p 2 ), (e x ®u 2 ), and (e x ®pu) which are effectively dealt with by 
considering the generalized characteristic function 

C=(e K+XQ ), (3.39) 

where K is a bilinear form in u and p 

K = ap 2 + pup + 7« 2 . (3.40) 

This procedure is equivalent to considering the Fokker-Planck equation for the joint distri- 
bution P(u,p,Q,t) as discussed in the previous subsection. The idea is to choose K, i.e., 
the parameters a, (3, and 7, in such a way that the unwanted correlations vanish yielding 
an equation for C. The conditions on K then yields the large deviation function fi directly. 
Embarking on the actual procedure below, we introduce the notation 

K' = K + AK, (3.41) 
Q' = Q + AQ, (3.42) 

where inserting ( 13.351) and ( 13.361) to order At 

AK =2ap(-(Y l + Y 2 )p -ku + F 1 + F 2 )At 

+(3{p 2 + u{-{Fi + T 2 )p -ku + F 1 + F 2 ))At 

+2-fupAt, (3.43) 

AQ =(F lP - T lP 2 )At + ^{F X F 2 + F 2 )(At) 2 , (3.44) 

Inserting in C = (exp(K' + \Q')) and expanding to O(At) we have 

C = (e K+XQ [l + AK + \AQ + ^(AK + XAQ) 2 }). (3.45) 
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Using the identity {F 2 exp(— F 2 /2A)) = A(exp(— F 2 /2A)) we can average over F\ and F2 
according to f!3.33|) and f l3.34|) inside the noise average defining C. We obtain after some 
algebra collecting terms to O(At) 

C = C + fiCAt + (e K+XQ {Ap 2 + Bpu + Du 2 ))At, (3.46) 

where the intermediate parameters A, B, D and \x in terms of at, (3, 7 and A are given by 

A = 4a 2 (riTx + T 2 T 2 ) + 2a(2Ar 1 T 1 - (T 1 + T 2 )) + f3 - XT 1 + \ 2 T 1 T 1 , (3.47) 

B = -2an - /3(T 1 + T 2 - 2XT l T l ) + 4a/3(r 1 T 1 + T 2 T 2 ) + 2 7 , (3.48) 

D = \p 2 - (3k, (3.49) 

H = 2a(T 1 T 1 + T 2 T 2 ) + XT 1 T 1 . (3.50) 

We note that the expression (13.461) involves correlations between exp(K + XQ) and p 2 , u 2 
and pu. However, since K is arbitrary we can obtain closure by choosing K, i.e., a, (3 and 
7, in such a manner that A = 0, B = 0, and D = 0. In the limit At — )■ (13.461) then reduces 
to the differential equation (I3.30p and /i locks on to the large deviation function 

In the present case of a bound Brownian particle the discussion is particularly simple. 
The condition D = immediately implies the two solutions /3 = and = 2k. However, 
since /x = for A = 0, the solution (3 = 2n must be discarded and we set (3 = 0. Likewise, 
7 is chosen so that B = 0. Finally, the condition A = yields a quadratic equation for a 
with admissible solution 

I\ + T 2 - 2Ar 1 T 1 - y /(T 1 + T 2 ) 2 + 2T 1 T 2 (1 - 2XT, + 2AT 2 - 2A 2 T!T 2 ) 

a(A) = 4( ri T 1 + r 2 T 2 ) ' (3 - 51) 



and we recover the case (11.5p for the free Brownian particle, i.e., 
MA) = \ 



r\ + r 2 - \ Y\ + r 2 + 21^2(1 - 2XT X + 2AT 2 - 2A 2 t\t 2 ) 



(3.52) 



IV. NUMERICAL SIMULATIONS 



Here we perform a numerical simulation of eqs. (12.1j) -( l2^2"j) . in order to sample the heat 
probability distribution function (PDF) P(Q,t) at long times and to verify that the distri- 
bution is independent of the spring constant k and in conformity with the large deviation 
function /i given by (jl.5p . Here and in the following the quantities will be expressed in 
dimensionless units. 
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FIG. 2. Large deviation function //(A) as a function of A, as given by eq. (|4.ip . for Ti = 1, r 2 = 2, 
T\ = 1, T 2 = 2. The shape is that of a half circle lying between the branch points A±, as given by 
T2D. 



Following Visco [35), see also 20|, |3jJ, /i(A) can be expressed in the form 

ri + r 2 



MA) 

where the branch points are given by 

1 



v / r 1 r 2 T 1 T 2 V(A + -A)(A- A- 



A- 



1_1± 1(1 1 V 1 ( r i + r ^) 2 

T x T 2 VV^i T 2 T 1 T 2 T 1 T 2 



(4.1) 



(4.2) 



note that A + > and A_ < 0. In Fig. [2] we have depicted the large deviation function /x(A) 
as a function of A. 

The large deviation function F(q), q = Q/t, characterizing the heat distribution, is 
determined parametrically from the large deviation function /j,(X) according to the Legendre 
transformation 



g = M(A) A* = A(g), 
F(g)=/i(A*)-AV(A*). 



We have, see also Visco 



35] 



F(q) = - [Tx + T 2 - q(X + + A_) - (A + - A_) v / r 1 r 2 T 1 T 2 + g2 



or inserting the branch points 



m = 5 



l l /(—- — 



1 1 y (r 1 + r 2 )^ ; , - 2 

rrnj + T^T 2 - VTiT2TiT2+q 



(4.3) 
(4.4) 

(4.5) 
(4.6) 
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FIG. 3. Heat PDF P{Q/t mSLX ) as a function of Q/t max for Ti = 1, T 2 = 2, T x = 1, T 2 = 2 and two 
different values of k : left panel i max = 10, right panel i max = 100. Full line: theoretical prediction 
as given by (|4.6p . Linepoints: PDF as obtained by simulating 10 5 independent trajectories. Inset: 
log- linear plot. 

Inspection of this equation shows that for small q we have a displaced Gaussian distribution; 
for large q we obtain exponential tails originating from the branch points \± in /i(A), i.e., 

F(q) ~ -\ + q for q > 0, (4.7) 
F(q) ~ -|A_||g| for q < 0. (4.8) 

In Fig. [3] we have depicted the distribution function P(Q/t) oc exp(tF(Q/t)), with F(Q/t) 
given by (I4.6p . as a function of Q/t on linear scales and log-linear scales (the inserts), for 
Ti = 1, T 2 = 2, T\ = 1, T 2 = 2, two different times t max = 10, 100, and two different values 
of the force constant k = 1, 10. We find good agreement between the simulations and the 
analytical results for the "central" part of the distribution. As expected, such an agreement 
improves as t max increases, being excellent for t max = 100. The tails cannot be sampled 
by the simulations, as they correspond to rare trajectories, that would require a very large 
simulation time to be observed. 

To further support our main finding, namely that the heat PDF is independent of the 
spring constant k, we calculated the first four moments of the distribution, over six orders 
on magnitudes of k, 10 -2 < k < 10 4 . The simulations were run for t max = 100, and 10 5 
independent trajectories were sampled. The results are reported in Fig. HI In the left panel 
we plot the relative change (Q m (K)) / (Q m (K = 0.01)), with m = 1 . . .4, and we find that 
the moments are practically constant over such a large range of values of k. Furthermore, 
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0.992 




1000 



FIG. 4. Analysis of the first four moments as obtained by numerical simulations with t max = 100, 
and 10 5 independent trajectories. Left panel relative change (Q m (K)) / (Q m {n = 0.01)) of the first 
four moments of the heat PDF as function of the spring constant k, wrt their value at k = 0.01. 
The moments are practically constant over a range of six orders of magnitude of k. Right panel: 
deviation of the first four moments from the expected value e m , as defined by (|4.9p . 



for each value of k, we calculate the deviation e m of such moments from the expected value 
which reads: 

(Qmim) ~ (Qex) 



(a m x ) 

where (Q™ um ) is the ra-th moment as obtained by the numerical simulations, and (Q™) is 
the corresponding exact value as obtained by equation (14. ip . The quantities e m are plotted 
in the right panel of fig. HI We find, that such deviations are negligible, basically due to 
numerical imprecision. 



A. Numerical investigation of the fourth-order potential case 

In the present subsection, we investigate the heat PDF of a particle coupled to the two 
heath baths at temperature T\ and T 2 , but moving in a quadratic potential 

V 4 (u) = a 2 u 2 + a 4 w 4 . (4.10) 

Thus in (12.21) the linear force is replaced by a term 2a 2 u + 4a^u 3 . We sample the heat PDF 
by considering 10 5 independent trajectories, with t max = 100, and choose different values for 
the parameters a 2 and in the potential (I4.10p . The results for the first four moments are 
reported in tabled and they provide a strong evidence that also in this case the heat PDF, 
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and so the large deviation function, is independent of the details of the underlying potential. 
As a bonus we also find that the first four moments are well described by the same large 
deviation function that we derived for the quadratic potential, which is independent of the 
potential details indeed, in the present case of the parameter a<i and a 4 appearing in (14. 10p . 



TABLE I. Deviation e m of the first four moments ( ()4.9p ) from the values predicted by the substrate- 
independent large deviation function, (|4.ip . The quantities 02 and 04 are the parameters of the 
fourth-order potential V4 as given by (|4. 10|) . 



a-2 


04 


ei 


£2 


£3 


£4 


-3 


1/2 


1.1 x 1(T 3 


1.9 x 1(T 3 


2.5 x icr 3 


3.0 x 10 -3 


-3/2 


1/12 


1.1 x 1(T 3 


1.6 x 1(T 3 


1.7 x 10~ 3 


1.6 x 10~ 3 


1 


1 


1.0 x 1(T 3 


1.3 x 1(T 3 


1.3 x 1(T 3 


1.0 x 10~ 3 



V. DISCUSSION AND CONCLUSION 



In this paper we have discussed a bound Brownian particle coupled to two distinct reser- 



voirs, generalizing a model proposed by Derrida and Brunet 31]. The issue was to determine 
whether the presence of a harmonic trap has an effect on the heat transport between the 
reservoirs and on the large deviation function characterizing the long time heat distribution 
function. By a variety of analytical arguments based on a Langevin equation evaluation 
of the two lowest cumulants and an evaluation of the large deviation function by a direct 
method due to Derrida and Brunet, supported by a numerical simulation, we have demon- 
strated that the presence of a harmonic trap has no effect on the heat distribution function 
which has the same form as in the unbound case. This result is maybe intuitively evident 
since a single particle, in contrast to an extensive system, does not have internal degrees of 
freedom. Furthermore, we provide numerical evidence, that the heat distribution function 
is unchanged if we consider a fourth-order potential, again supporting our finding that such 
a distribution is independent of the underlying potential. 



It also follows that the Gallavotti-Cohen fluctuation theorem 19] in ( 11.21) is unchanged 
by the presence of the spring. The fluctuation theorem is associated with the entropy pro- 
duction Q1/T1 and Q2/T2 at the heat sources whereas the presence of the spring represents 
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a deterministic constraint not associated with entropy production 
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